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ABSTRACT 

The equilibria formed by the self-gravitating, collisionless collapse of simple initial conditions have 
been investigated for decades. We present the results of our attempts to describe the equilibria formed 
in A^-body simulations using thermodynamically-motivated models. Previous work has suggested that 
it is possible to define distribution functions for such systems that describe maximum entropy states. 
These distribution functions are used to create radial density and velocity distributions for comparison 
to those from simulations. A wide variety of A-body code conditions are used to reduce the chance 
that results are biased by numerical issues. We find that a subset of initial conditions studied lead to 
equilibria that can be accurately described by these models, and that direct calculation of the entropy 
shows maximum values being achieved. 

Subject headings: galaxies:structure — galaxies:kinematics and dynamics 


1. INTRODUCTION 

Over the past several decades, numerous investigations 
of collisionless, self-gravitating systems have been under¬ 
taken. Erom early focus o n the formation an d evolution 
of elliptical galaxies fe.g.. Ivan Albadalll982h . to cosmo- 
logic al simulations of dark matter structure formation 
(e.g.-lNavarro Frenk. fc WhitdllQQ^ iMoore et al.l 119991 : 
iSpringel et al.1120051 1. these works have taken advantage 
of ever-increasing levels of computing power. Simula¬ 
tions with higher and higher resolutions (numbers of par¬ 
ticles per unit volume) are constantly being performed. 
Our goal is not to attempt to replicate these state-of- 
the-art simulations, but rather to use more modest sim¬ 
ulations to investigate some basic questions. The sys¬ 
tems that we simulate are not direct analogues of puta¬ 
tive dark matter halos nor elliptical galaxies, but they 
do share the fundamental physical conditions of self¬ 
gravitation and collisionless evolution. As many previ¬ 
ous works have shown evidenc e of “universal” behaviors 
[such as radial deiisity ( e.g., iNavarro. Frenk. fc White! 
Il997t iNavarro et ^120041) and powe r-law “phase space” 
p{r)/a^{r) (iTavlor fc: Navarroll2001h profiles], an inter¬ 
esting possibility is that a basic physical mechanism may 
underlie the for mation of these self-gravitating e quilibria. 

Earlier work (|Barnes fc WilliamsI 20 III 120121) presents 
descriptions of distribution functions of collisionless, 
self-gravitating systems that represent maximum en- 
tropv states . The se results follow the seminal work of 
iLvnden-Belll 1)19671 ) , where it is argued that a fourth sta¬ 
tistical family is appropriate for describing the phase- 
space evolution of these kinds of systems. In a nut¬ 
shell, the familiar Maxwell-Boltzmann statistics describe 
systems in which particles are distinguishable and do 
not obey a phase-space exclusion principle - multi¬ 
ple particles can occupy a very small region of phase- 
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space. On the other hand, the Lynden-Bell statisti¬ 
cal family is appropriate for systems of distinguishable 
particles that follow an exclusion principle ~ a classi¬ 
cal version of Fermi-Dirac statistics. By relaxing the 
requirement of large pha se-space occupation numbers, 
iBarnes fc WilliamsI (|2012[ ) show that finite-mass, max¬ 
imum entropy states exist for both Lynden-Bell and 
Maxwell-Boltzmann statistical families. The major goal 
of this work is to test whether or not any simulated 
system will relax to such a state. It is certain that 
these models will fail for sufficiently strong collapses, 
as the assumed velocity isotropy underlying the mod¬ 
els has to disappear as the ra dial orbit instability be¬ 
gins to become important (e.o.. iMerritt fc Aguilad(l985t 
IBarnes. Lanzel. fc Williamsll2009t) . As such, we also aim 

to identify ranges of conditions that can result in maxi¬ 
mum entropy states. A final goal is to monitor the be¬ 
havior of entropy in simulations to see if a maximum 
value is reached. 

We have created suites of simulations to test the use¬ 
fulness of the Lynden-Bell and Maxwell-Boltzmann fam- 
ilies of models. The publicly-available GADGET code 
(|Springell I2005D has been employed to evolve simula¬ 
tions of collisonless systems comprised of A = 10® and 
A = 10® particles. These simulations approximate colli¬ 
sionless conditions by utilizing softened interactions. As 
a point of comparison, we have also analyzed collisional 
systems with A = 2^^ Ri 10® particles using a version of 
NBOD Y-6, enhanced with a Gra phics Processing Unit 
(GPU) (|Nitadori fc Aarse'thll2012h . This code uses direct 
Newtonian particle-particle interactions, but the large 
number of particles guarantees that two-body relaxation 
processes occur over times^ cales thousands of times longer 
()Binnev fc Tremainelll987l Ch. 4) than any gravitational 
potential (or “violent”) relaxation processes, which occur 
over a few initial crossing times T. 

For any given simulation, we analyze spherically sym¬ 
metric density and velocity distributions by counting par- 
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tides in spherical bins centered on the system center-of- 
mass. With estimates of the various uncertainties, these 
radial profiles are then used as the data to be matched by 
the Lynden-Bell and Maxwell-Boltzmann models. Chi- 
squared minimizations are used to indicate the appropri¬ 
ateness of each model. We will not insist that models 
with low chi-squared values are the only descriptions of 
these simulations, merely that such a model is consistent 
with the data. For simulations that are well-described 
by the thermodynamic models, we also investigate the 
behavior of entropy during its evolution. 

We begin by describing simulation initial conditions, 
evolution code details, and analysis techniques in sec¬ 
tion Two methods for describing entropy behavior in 
these simulations are presented in Section [3] Section |4] 
contains the findings inferred from minimizations for se¬ 
lected simulations, while section [5] outlines the entropy 
behaviors seen in the simulations. We summarize in Sec¬ 
tion [51 

2. SIMULATION DETAILS 
2.1. Initial Conditions 

For all simulations discussed here, particles are as¬ 
sumed to be identical and system mass, system radius 
and Newton’s gravitational constant are set to unity 
(msys = R = G = I). Each sy s tem i s composed of 
N particles. Following Ivan Albadal (|I982ll . particle posi¬ 
tions are chosen according to two different schemes. 

In “single” simulations, initial particle locations are 
randomly chosen within the system according to a spec¬ 
ified density distribution ~ cuspy (p oc 1/r) or Gaussian 
(p oc exp—r^). A simple rejection scheme is used to 
generate the distributions. The system center-of-mass 
must coincide with the center of the spherical boundary 
to better than I/-\/]V to be acceptable. 

For “clumpy” simulations, centers-of-mass locations 
are chosen for small clumps of particles according to the 
above density distributions, and then particles are uni¬ 
formly distributed within a clump. The numbers of par¬ 
ticles in each clump are chosen from a Salpeter distribu¬ 
tion - the probability of generating a clump with Ndump 
particles is proportional to where a = 2 jZ for 

this work. As we are not investigating any specific phys¬ 
ical situation with these simulations, the details of this 
distribution (for example, variations to the adopted a) 
have not been scrutinized. The Salpeter form has been 
adopted based on its simplicity. We demand that the sum 
of the volumes of all the clumps equals that of a sphere 
with radius R = 1. Individual clump radii are chosen 
proportional to the fraction of the total mass they con¬ 
tain, rdump c)c (A"ciump/fV)^^^- Clumps can overlap one 
another, leading these initial conditions to have regions 
of higher-than-average density as well as nearly empty 
regions where clumps fail to overlap. 

Particle velocities are given random orientations, guar¬ 
anteeing initial velocity isotropy. Speeds are chosen by 
adopting an initial virial ratio Qo = 2Aro/|Wo| that links 
a system’s initial kinetic energy Kq to its initial potential 
energy Wq. Once particle positions have been selected, 
the virial ratio is used to define a scale speed. For single 
simulations, this is the speed given to every particle in 
the system. Clumpy simulations distribute speeds in a 
more complicated manner. We define hot-clumpy sys¬ 


tems to be ones in which the clump centers-of-mass have 
zero initial velocities - particles in clumps are given ran¬ 
dom velocity directions with equal speeds. Cold-clumpy 
systems are composed of clumps in which all of the par¬ 
ticles move with the clump center-of-mass velocity. The 
centers-of-mass velocities are randomly oriented but have 
the same magnitudes. As an intermediate case, warm- 
clumpy simulations split the kinetic energy equally be¬ 
tween individual particle motion and clump center-of- 
mass motion. Independent of the specifics of the setup, 
the systems are not initially in mechanical equilibrium, 
even though they may be in virial equilibrium (if Qq = 1 
is adopted). 

2.2. Evolution Code Details 

As mentioned in the introduction, this work utilizes 
two very different codes for evolving initial conditions. 
Our aim is to be able to identify any numerical effects 
due to particle number, softening parameters, and/or 
code specifics. The GPU-enhanced NBODY-6 code has 
an architecture designed for investigating globular clus¬ 
ter dynamics. GADGET has been designed to perform 
cosmological simulations of structure formation. Find¬ 
ing agreement between our predictions and the results of 
both types of simulations will strengthen the assertion 
that our analytical picture is relevant to collisionless sys¬ 
tems and is not biased by numerical issues. 

GADGET is a versatile tree-code that incorporates 
softened forces. In general, we have adopted the standard 
parameters for GADGET. However, we do not evolve in 
an expanding universe, and we adopt different softening 
lengths, depending on the situation. For N = 10® simu¬ 
lations, we adopt a softening length e = 10“^, about 100 
times smalle r than the “optima l” softening length value 
described in iPower et al\ (120031 1. Test runs with soften¬ 
ing lengths between our adopted value and the optimal 
value have resulted in only minor differences in density 
profile shape. Smaller values lead to integration times 
that we judged to be unacceptably long, so our value 
is as small as possible while keeping the wall-clock time 
manageable. For N = 10^ simulations, we adopt the op¬ 
timal softening length e = 4 x 10“®. Again, testing with 
smaller softening lengths indicates that density profiles 
are largely unaffected, but integration times significantly 
lengthen. 

We have performed three types of GADGET simula¬ 
tions. For N = 10®, we have varied Qo between 0.1 
and 1.0 for both types of initial density distributions and 
for the three different velocity assignments when start¬ 
ing with clumpy initial conditions. We take these as 
the standard set of simulations that provide zeroth-order 
tests of our models. As the initial conditions are based 
on random distributions of particles, we have also per¬ 
formed ensemble simulations of a subset of these initial 
conditions. Five independent realizations of initial con¬ 
ditions with 0.7 < Qo < 1.0 and both initial density dis¬ 
tributions have been evolved. For clumpy simulations, 
only ensembles with hot conditions were created. Sec¬ 
tion 0] discusses the justihcation for the ranges and spe¬ 
cific values used. Averages of these five realizations have 
been used to validate the results of the standard simu¬ 
lations. The third type of GADGET simulation involves 
fV = 10® particles. Both initial density distributions with 
0.7 < Qo < 1-0 have been investigated for single and 
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clumpy systems. Like the ensemble simulations, these 
additional results provide an estimate of the robustness 
of the results based on the N = 10® simulations. 

The GPU-NBODY-6 simulations evolve single and 
clumpy systems with N > 10®, 0.7 < Qo < 1.0, and 
both initial density profiles. This code does not utilize 
softened forces, so it is a collisional code. However, the 
large number of particles provides a reasonable basis for 
the assumption that two-body effects {e.g., ejection) have 
relatively minor impact on the simulations. For exam¬ 
ple, no particles gain escape speeds during any of our 
simulations due to two-body collisions. Unfortunately, 
simulations with larger particle numbers could not be 
completed with the current hardware available to the au¬ 
thors. 


2.3. Analysis 

Independent of the evolution code, each simulation ex¬ 
tends for at least 20 initial crossing times. We have ob¬ 
served that this is generally a sufficient period for a sys¬ 
tem to reach a mechanical and virial equilibrium state. 
For clumpy initial conditions, individual clumps have 
“dissolved” by 10 initial crossing times, at the latest. By 
construction, clumps in the hot-clumpy simulations be¬ 
gin to disperse almost immediately. Any simulation with 
a longer evolution will be noted in what follows. We in¬ 
fer mechanical equilibrium by observing that the radius 
of the inner-most 90% of the particles stops varying and 
that the average velocities of the particles are zero over 
the radial range of the inner-most 90% of particles. Sim¬ 
ulations with Qo ^ 0.3 typically reach these conditions 
in less than 20 initial crossing times. Systems are decom¬ 
posed into spherical, 1-percent mass shells at each output 
timestep. Particles in these shells provide density and 
velocity statistics (averages, rms values, and variances). 
Each shell is broken into three sub-shells which provide 
ranges for the density and velocity values that are used 
to estimate uncertainties (maximum minus minimum). 
Systemic axis ratios, phase-space occupation values, and 
entropy production rates (see section El) are also calcu¬ 
lated at each output timestep. 

With the radial density p and rms speed tirms distribu¬ 
tions formed, we next compare to several model distri¬ 
butions. Our focus is on the comparison to the Lynden- 
Bell and Maxwell-Boltzmann models, but we also include 
two common analytical models, Plummer and de Vau- 
couleurs. The Plum mer models we use have a density dis¬ 
tribu tion given by (iPlummed 119111 : iBinnev &: Tremainel 

[HU, 



where po is a scaling density and tq is a scaling radius. 
This density profile is relatively constant in the core of 
the system and declines rapidly ( oc r~®) near the outer 
edge. The de Vaucouleurs profile (Ide Vaucouleursl[l948[) 
has been used to fit the light profiles of elliptical galax¬ 
ies, and is a simple example of a broken power-law distri¬ 
bution (a cousin to commonly- discussed, cosmologically- 
motivated profiles such as, iNayarro. Frenk. fc Whit j 
Il996t iMoore et al.l 1199911 . It is well-established that 
the type of simulations discussed here result in struc¬ 
tures with outer density behayior that matches the de 


Vaucouleurs profi le when Qo ^ 0.2 (jyan Albadal Il982[ 
iSylos Labinill20l3l . so it proyides a good benchmark for 
our standard simulations. The de Vaucouleurs density 
profile is giyen by. 



where i5 = i, and po and tq are again a scaling density 
and radius, respectively. A de Vaucouleurs density profile 
has a central cusp, in contrast to the central density core 
behavior of the Plummer model. 

As our models and simulated systems all have finite 
masses, we demand that their density and Vrms values 
match at the half-mass radius. With the connection be¬ 
tween model and simulated data values fixed, we use the 
reduced chi-squared statistic as the figure of merit for 
our fits. 


1 (M, - A)^ 

A^data 


(3) 


where Vdata = 100 for our 1% mass shells, Di is a mass 
shell density or Uj-ms value from a simulation. Mi is a 
model value corresponding to the same radial location, 
and Ai is an uncertainty estimate for the simulation value 
(as described in the first paragraph of this section). One 
should expect, if the uncertainty estimates are appropri¬ 
ate, that a good model fit to the data produces Xr ~ 1- 

The Plummer and de Vaucouleurs density models have 
no free parameters, so their Xr values are determined 
upon matching to simulation values at the half-mass ra¬ 
dius. For Lynden-Bell and Maxwell-Boltzmann models, 
a single parameter (z^lb or izmb) determines the shape 
of the density profile, and hence the match to the sim¬ 
ulation. Density profiles are determined by iteratively 
solving the Poisson equat ion in straightforward fashion 
(|Binnev fc Tremainel 19871 Sec. 4.4.2). The best-fit value 
of v is determine d using an amoeba Xr minimization 
([Press et al.l 119941 ). We have also performed Markov 
Chain Monte Carlo minimizations to corroborate the 
amoeba results. 

To determine model Vmis distributions, we solve the 
Jeans equation for a given density prohle. This ap¬ 
proach demands a choice be made regarding the radial 
behavior of the velocity anisotropy f3(r). We utilize two 
Vrms profiles: one assuming velocity isotropy f](r) = 0 
and one adopting /3(r) from the simulated system. The 
Lynden-Bell and Maxwell-Boltzmann models are derived 
assuming velocity isotropy, so their density profiles are 
consistent only with the /3{r) = 0 ?;rms profiles. In the 
absence of a distribution function that incorporates the 
mild tangential velocity anisotropy present in our sys¬ 
tems, we assume that the density derived from such a 
function should be well-approximated by the density re¬ 
sulting from the isotropic version of the distribution func¬ 
tion. While not an exact description of the situation in 
our work, this assumption is consistent with the behavior 
of the aniostropic Plummer model described in IMerritd 
(|198,^ . We choose the /3(r) to be used in the Jeans equa¬ 
tion to be a smoothed version of the velocity anisotropy 
present in the simulation. Specifically, a fourth-order 
polynomial fit to a simulation anisotropy profile is cre¬ 
ated. The parameters describing the Lynden-Bell and 
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Maxwell-Boltzmann models are not allowed to vary dur¬ 
ing comparison to the tirms distribution, so the velocity 
calculation for all models is straightforward once the 
model and simulated Urms values are matched at the half¬ 
mass radius. 


3. ENTROPY BEHAVIOR 
3.1. Microscopic Picture 

The basis of the iBarnes fc WillianisI (|201lL l2012f) 
work is the phase- space counting approach outlined in 
iLvnden-Belll (|1967ll . Here, we present a brief summary 
of the chief ideas necessary for defining entropy from 
this microscopic viewpoint. Phase-space is imagined to 
be sub-divided into two lattices; an array of nearly in¬ 
finitesimal micro-cells (each with volume w) that are ar¬ 
ranged into collections of macro-cells. Macro-cells con¬ 
tain V micro-cells, and this value serves as the control 
parameter for the models. Micro-cell occupation defines 
a fine-grained distribution function, while macro-cell oc¬ 
cupation defines a coarse-grained distribution function 
that can be realized through simulations. The occupa¬ 
tion of micro-cells determines the statistical properties 
of the system. Maxwell-Boltzmann statistics arise when 
multiple distinguishable particles can occupy a micro¬ 
cell. Lynden-Bell statistics describe situations in which 
a classical exclusion principle disallows multiple particles 
occupying a single micro-cell. For collisionless systems 
like the ones we investigate here, the fine-grained distri¬ 
bution function is a constant of motion. As such, one 
expects the Lynden-Bell statistical family to be most ap¬ 
propriate because the fine-grained distribution function 
values cannot be increased through multiple occupancy. 

With these ideas, one can count the number of en¬ 
ergy states a vailable to a s ystem , and hence, define the 
entropy. The iLvnden-Belll (Il967ll work relies on the Stir- 
ling; approximation to sim plify entropy calculations, but 
IBarnes fc Williainsl (|2012ll relax this assumption and al¬ 
low macro-cell occupation numbers to be small enough 
that the Stirling approximation fails. The end result is 
that the Lynden-Bell entropy c an be given as (Equation 
13 in iBarnes fc Williamsll201^ . 


5'lb = •S'lb.o— 

M 

^B ^ [{ni + 1/2) In {rii -|- 1)-|- 

i=l 

iy — rii -\- 1/2) In {ly — rii + 1) + Xo,ni + i(4) 

where S'lb.o = fee [V In TV — V -|- M(z/ In v — hi 27r)], N is 
the number of phase-space elements/particles, and M is 
the total number of macro-cells. The A function arises 
from the approximation. 


Inx! = (a;-I-^) In (a;-I-1) - a;-I--I-Ao,a; (5) 


where 




(^^ + 2a: + ii) 

(x2 + f|a;+^)- 


( 6 ) 


A similar expression can be found for the 
Maxwell-Boltzmann ent ropy (see Equation A2 in 
IBarnes fc Williamsll2012f) . Adopting these modifications 
leads to the possibility of finite-mass and energy systems 


that belong to the Maxwell-B oltzmann statistical family, 
in contrast to the fin dings of iLvnden-Bell (1^671). Like¬ 
wise, the discussion in iBinnev fc Tremaind ( 1987L §4.7.1) 
regarding the impossibility of maximizing entropy be¬ 
comes invalid, as the simple / In / term in the entropy 
calculation is now modified. Our Lynden-Bell expression 
does not change the overall charac ter of the as s ociate d 
distribution function presented in ILvnden-Belll (Il967ti : 
it remains finite-mass and energy and closely resembles 
a Fermi-Dirac distribution. 

We have calculated Slb using the results of GADGET 
N = 10® simulations. At every timestep, positions and 
velocities are used to assign each particle to a macro-cell, 
giving Hi. A fixed value of z/lb = 10“* has been chosen 
for these calculations, as that is always greater than the 
maximum rii value in these simulations. Tests varying 
z/lb show that the “zero point” of S'lb is affected much 
more strongly than its time-dependent behavior . Results 
of these calculations are discussed in Section [5] 


3.2. Macroscopic Picture 

As a complement to the microscopic approach, we fol¬ 
low the discussion of thermal non- equilibrium situations 
given bv ide Groot fc MazuH (jl984ll . The behavior of self- 
gravitating systems composed of a large number of mas¬ 
sive particles can be described using equations that rep¬ 
resent macroscopic conservation laws. Expressions of the 
conservation of energy can be manipulated and combined 
with the first law of thermodynamics to provide insight 
into the behavior of entropy. In particular, it is useful to 
write the entropy production rate of a system as. 


dt 



(7) 


where tr is the entropy production per unit volume per 
unit time and the integral is taken over the system vol¬ 
ume. The specific entropy production rate in this picture 
is, 

CT =VTk - ^H : Vvo, (8) 


where Tk is the kinetic temperature, q is the heat con¬ 
duction flux, H is the anisotropic pressure tensor, and 
Vq is the mean velocity. As usual, the kinetic tempera¬ 
ture is a measurement of the random kinetic energy in 
a small region; Tk = (m/3kB)(Vp}, where ks is Boltz¬ 
mann’s constant and Vp is the magnitude of the peculiar 
velocity. The heat conduction flux q = §{VpVp) repre¬ 
sents the peculiar kinetic energy that is transported by 
peculiar velocities in the system. Int erested readers may 
hnd d etails of the calculation of a in Ide Groot fc MaznrI 
dun §3.3). 

In order to attempt to follow the behavior of entropy 
from this macroscopic viewpoint during a simulation, we 
form the necessary macroscopic quantities by averaging 
over spatial volumes. Simulated systems are divided into 
spherical volume elements, and values for temperature 
and pressure tensor components are found using parti¬ 
cles within the volume element. While this procedure is 
straightforward, the average number of particles per ele¬ 
ment can be rather small, even for modest spherical grid 
resolutions. For example, in a simulation with N = 10® 
particles broken into a spherical grid with 10 radial, 8 
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polar, and 8 azimuthal bins, each will contain on the or¬ 
der of 100 particles. The need for gradients in quantities 
places some constraint on how coarse the spherical grid 
can be made. We use simulations with N = 10® parti¬ 
cles and the aforementioned grid resolution to begin to 
investigate this macroscopic picture of entropy behavior 
during an evolution. Results of these calculations are 
discussed in Section [5] 

4. DENSITY AND VELOCITY PROFILE FITTING 

We take the GADGET simulations with N = 10® as 
the standards for the majority of our results. These simu¬ 
lations provide a wide-range of evolutions from which we 
draw broad inferences. The other simulations (ensemble 
N = 10® GADGET, N = 10® GADGET, GPU-NBODY- 
6) are focused on testing relationships and behaviors sug¬ 
gested by the standard set. 

4.1. Individual N = 10® GADGET Simulations 

Lynden-Bell and Maxwell-Boltzmann models can pro¬ 
vide excellent fits to the density and velocity profiles of 
the final equilibrium states of single and hot-clumpy sys¬ 
tems when Qq > 0.6. However, independent of the type 
of particle distribution (single or clumpy) and initial den¬ 
sity profile, systems with Qq < 0.6 do not evolve to states 
like those predicted by our maximum entropy argument. 
The density and velocity distributions of these simu¬ 
lated equilibria show significant deviations from thermo¬ 
dynamic model expectations. For systems with small 
enough Qq and/or cold-clumpy initial conditions, the fi¬ 
nal density distributions are best-described by a de Vau- 
couleurs profile, due to its cuspy nature. However, the 
specifics of the central density cusp do not always agree 
with the 6 = 1/2 value of the de Vaucouleurs profile. We 
have effectively let 6 vary in a few cases, finding that 
1/4 ^ 6 ^ 3 / 4 . The thermodynamic models investi¬ 
gated here simply cannot reproduce the central density 
cusp. Interestingly, Plummer models provide good de¬ 
scriptions of the density profiles of warm-clumpy systems 
when Qo > 0.6. 

Now, we turn to the simulations where there is bet¬ 
ter agreement between the thermodynamic models and 
the simulations. As illustrations of the quality of these 
fits. Figures [T] and [2] contain plots of logarithmic den¬ 
sity profiles for single and hot-clumpy equilibrium sys¬ 
tems evolved from both initial density profiles. In each 
panel, the LB and MB models are superimposed (com¬ 
parison Plummer and de Vaucouleurs profiles are also 
shown in the upper-left-hand panel). To highlight the 
range of density profile behaviors. Figure [T] contains re¬ 
sults of simulations with Qq = I.O, while Figure [2] shows 
results from evolutions with Qo = 0.7. Analogous plots 
for simulations with Qq = 0.8 and Qq = 0.9 (not shown) 
reveal very similar model behaviors. For single systems, 
LB models produce fits superior to those from MB mod¬ 
els in every case. For hot-clumpy systems, MB models 
perform better than LB models for Qq > 0.9, with the 
models reversing positions for Qq < 0.8. Overall, LB 
models appear to do a better job of describing the den¬ 
sity distributions of the simulated equilibria. 

Corresponding plots for the Urms profiles are given in 
Figures |3] and Si In these figures, the profiles created 
assuming isotropic and anisotropic velocity distributions 


show quite different behaviors. Finding Vrmsif) from the 
Jeans equation under the assumption that /3(r) = 0 re¬ 
sults in profiles that have flat cores. However, single sys¬ 
tems generically have Urms profiles with non-zero slopes 
near their centers. Hot-clumpy systems are qualitatively 
similar to the isotropic model predictions, but are not 
terribly well described by the models. The set of thin 
lines below the data points illustrate the Urms.r, "Crms.e) 
and Vrms,(j> behavior of the simulation. The differences 
between the line shapes indicates mild tangential veloc¬ 
ity anisotropy. Allowing this /3(r) profile as input to the 
Jeans equation results in Umis profiles that dramatically 
increase the quality of the fits. We note that the LB 
models tend to provide better fits to the simulated v^ms 
profiles than those produced by the MB models. 

4.2. Other Simulations 

Fits to the standard simulations suggest that LB mod¬ 
els are generally better than MB models. We now begin 
to test the robustness of this observation using “aver¬ 
age” systems formed by ensembles of simulations with 
different realizations of the same initial conditions using 
N = 10® particles. As with the standard simulations, 
GADGET has been used to evolve the initial conditions. 
Figures [5] and [ni are analogous to Figures [T] and [2] The 
most significant changes one notices is that the simula¬ 
tion profiles are smoother and have smaller error bars. 
For nearly every simulation, LB model fits to the den¬ 
sity distributions are superior to those provided by the 
MB model (for Qo = 1, the two models can produce 
comparable fits). The same is true for velocity profiles. 
Figures [7] and IH] show the improvement that LB models 
provide over MB models. As with the standard simu¬ 
lations, the inclusion of velocity anisotropy significantly 
improves agreement between the model curves and the 
simulation results. 

The equilibrium structures in simulations with N = 
10® particles (returning to only one realization per ini¬ 
tial condition) have also been analyzed in a similar man¬ 
ner. Figures 1^ and ITOl show how the LB and MB density 
profiles compare to the simulation results. The single 
system profiles (particularly for the cuspy initial profile) 
display small-scale variations that are larger than the 
uncertainty estimates we have made. The systems are 
in virial equilibrium, and the near-zero average velocity 
values for all components indicate mechanical equilib¬ 
rium as well. Additional evolutions of these initial condi¬ 
tions with different (smaller) softening lengths have pro¬ 
duced very similar outcomes. Extending the evolutions 
to longer times does reduce the variations somewhat, and 
we have used the results of our longest evolutions (out to 
30r) to create the relevant figures. Somewhat surpris¬ 
ingly, we do not see similar structures in the profiles of 
the hot-clumpy simulations, which are very smooth. As 
these new features suggest that our naive uncertainty es¬ 
timates are questionable, we do not place much stock in 
the actual values of Xr that we have found. However, we 
argue that since the LB and MB models are both being 
compared to the same data with the same uncertainties, 
their relative Xr values distinguish between which is the 
more appropriate model. 

Unlike the averaged simulations, the appropriateness of 
the LB model is not obvious here. In most cases, the LB 
and MB models provide comparable fits to the simulation 
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density profiles; only for the hot-clumpy simulations with 
Gaussian initial density profiles is the LB model clearly 
preferred. Likewise, the differences between the LB and 
MB velocity profiles are more subdued in Figures [TT] and 
1121 compared to previous versions. It also appears that 
there are more significant differences between the data 
and the LB model curves in these figures. 

Overall, the GADGET-based simulations support the 
idea that the LB model does a good job describing the 
density and velocity profiles of the equilibria of colli¬ 
sionless systems that have undergone mild gravitational 
potential relaxation. Our final test for this idea is to 
use the different evolution code GPU-NBODY-6. Unlike 
GADGET, GPU-NBODY-6 does not incorporate grav¬ 
itational softening and instead treats particles as true 
point masses. While two-body encounters do occur, the 
large particle number (TV > 10^) minimizes the global im¬ 
pact of particle-particle interactions, leaving a basically 
collisionless evolution. The results of these simulations 
are similar to those from the standard simulations. In 
general, LB and MB models provide comparable descrip¬ 
tions of the density profiles of systems with Qq > 0.9, 
but LB models are superior for Qq < 0.8 (see Figures [13] 
anddH). Figures [T5] and [T51 show that LB models tend 
to provide better fits to Vrms profiles produced by these 
simulations. 

5. ENTROPY PRODUCTION 

As mentioned in section ICT the microscopic calcula¬ 
tion of the entropy has been carried out for the TV = 10® 
GADGET simulations. We have created two different 
macro-cell grids; one with 10 macro-cell divisions per 
phase-space dimension {M = 10®), and one with 15 di¬ 
visions {M K, 10^). We show a representative pair of 
ShB{t) curves in Figure [T7h.. The different curves corre¬ 
spond to the grid choices indicated. These are derived 
from a single, cuspy Qo = 0.7 system, but similar behav¬ 
ior is seen in simulations with other Qo, Gaussian density 
profiles, and clumpy particle distributions. Most impor¬ 
tantly, it is clear that the entropy rises from an initial 
value and rapidly approaches a maximum, steady-state 
value. The most rapid increase (a roughly 5% change) in 
entropy occurs during the first couple of initial crossing 
times. It seems natural that the most rapid growth in 
S'lb occurs during the period when the strongest gravi¬ 
tational potential relaxation (biggest variations in poten¬ 
tial) occur. For comparison, the variation in the virial 
ratio Q = 2Ar/|lU| as a function of time in the same 
simulation is shown in Figure [T7b . 

We also note that the finer macro-cell grid produces 
smaller variations in the value of S'lb, a trend that also 
occurs as Qo 1- With the increasing smoothness of the 
curves with higher Qo, one can also discern that there 
is slower growth of Slb occurring over tens of cross¬ 
ing times. We speculate that this slower growth is at¬ 
tributable to the phase mixing that continues after the 
initial potential variations decrease. 

While the microscopic entropy picture dovetails nicely 
with the overall scenario of entropy production in colli¬ 
sionless systems, the macroscopic picture results are not 
so clear. Figure [THI illustrates the behavior of Equation^ 
for the single, cuspy Qo = 0.7 simulation. Given the mi¬ 
croscopic results, one would expect a rather tall, positive 
spike to appear in the early part of the evolution, fol¬ 


lowed by a decline towards zero. Instead, the very noisy 
curve seems to oscillate about zero. A boxcar smoothing 
filter applied to the raw values makes the oscillation more 
plain. Again, this same behavior is seen across the vari¬ 
ous TV = 10® GADGET simulations. Given the reliance 
on derivative information required by this approach, the 
noise present in the values is not surprising. It is possible 
that our particle number and grid resolutions are simply 
not high enough to capture the true behavior of the en¬ 
tropy creation term. Preliminary testing with higher grid 
resolution did not produce appreciably different results. 
However, there could also be a fault in the assumption 
of local thermodynamic equilibrium that underlies the 
derivation of Equation [H] 

6 . GONCLUSIONS 

We have evolved sets of TV-body initial conditions to 
determine if the thermodynamically-motivated Lynden- 
Bell or Maxwell-Boltzmann distribution functions can 
describe the equilibrium distributions of particle loca¬ 
tions and velocities. Different evolution codes and nu¬ 
merical parameters have been adopted to reduce the like¬ 
lihood of spurious findings. These simulation results have 
also been used to investigate the entropy behavior of 
these systems. 

Initially “cold” systems evolve in expected fashion, 
forming equilibria that show cuspy central density pro¬ 
files and outer density profiles that match the de Vau- 
couleurs form, not the thermodynamic models. How¬ 
ever, a subset of our simulations can be well-described by 
distribution functions that maximize entropy. Initially 
“hot” systems form equilibria with density profiles that 
are very similar to those produced by the Lynden-Bell 
distribution function. Unfortunately, the lack of velocity 
isotropy in these simulations seems to preclude the ability 
of the LB distribution functions to predict v^ms profiles. 
However, accounting for the mild tangential anisotropy 
produces extremely good descriptions of simulation ve¬ 
locity results. Previous simulations involving mild ve¬ 
locity anisotropy indicate ver y similar global behavio r 
to fully isotropic systems le.o.. iMerritt fc AgnilarlflO^ . 
We argue that the isotropic behavior is of central impor¬ 
tance, with velocity anisotropy determining high er-order 
corre ctions to density and Wrms profiles {e.g., iMerrittI 
1198511 . To be clear, introducing velocity anisotropy out¬ 
side the distribution function, as we have done, means 
that such models are not self-consistent. The anisotropic 
frms profiles are not predictions of the maximum entropy 
argument, which assumes only constant system mass and 
energy. Anisotropic models would require entropy max¬ 
imization in cluding at least one o ther constraint, along 
the lines of (jTrenti fc Bertinir2005[l . Given these caveats, 
our results suggest that there are conditions under which 
a collisionless self-gravitating system will evolve to a 
maximum entropy state. This conclusion is strengthened 
by the fact that direct calculation of the entropy (using 
a phase-space occupation approach) also shows a maxi¬ 
mum value being attained. An alternative construction 
of the entropy behavior is inconclusive, presumably due 
to resolution effects. 

There remain some unresolved issues with the idea of 
these equilibria representing maximum entropy states. 
One is that maximized entropy is normally taken to im¬ 
ply thermodynamic equilibrium, but these simulated sys- 
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terns clearly have kinetic temperature gradients. Is it 
possible that in self-gravitating systems, these two condi¬ 
tions are not equivalent? Another set of questions revolve 
around the role of the v parameter. This value represents 
the number of micro-cells that occupy any macro-cell, 
but it does not have an ab initio value. In the direct 
calculation, it must be larger than the largest macro¬ 
cell occupation number rii in order for the entropy to be 
well-defined. In fitting density profiles, we have placed no 
such restriction on its value. For the simulations focused 
on in this work, we have found 5 ^ 2000. In general, 

higher Qo values couple to lower v values. On the other 
hand, the slight positive concavity seen in the outer den¬ 


sity prohles of Qo = 0.7 simulations is reproduced by 
the models only when v > 1000. Unsurprisingly, there 
is no simple answer underlying the evolution of collision¬ 
less systems, but our results suggest there are some new 
questions to ask. 
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Fig. 1.— Logarithmic density profiles for individual N = 10^ GADGET simulations with Qq = 1.0. The radial coordinate is scaled 
by the half-mass radius (indicated by the vertical dashed line) in all profiles. In each panel, the initial conditions are specified and the 
errorbars indicate the data values from the simulations. The curves show the behaviors of the various models, specified in the legend. As 
they provide poor descriptions of the simulations, the Plummer and de Vaucouleurs models are only included in the single cuspy panel for 
reference. Lynden-Bell (LB) model fits produce density Xr values smaller or comparable to those with Maxwell-Boltzmann (MB) models 
(except for the hot-clumpy Gaussian simulation): single cuspy - = 0.635, Xmb “ 0-635; single Gaussian - X^b ~ 0.811, Xmb “ 2.162; 

hot-clumpy cuspy - Xlb “ 1-207, Xmb ~ 0.554; hot-clumpy Gaussian - Xlb “ 5.093, Xmb “ 0.581. 
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Fig. 2.— Logarithmic density profiles for individual N = 10^ GADGET simulations with Qo = 0.7. The panels and linestyles are the 
same as in Figure^ In these simulations, the density profiles share a slight upward concavity in their outer profile. In general, the LB 
model fits are able to match this data behavior better than the MB models: single cuspy — ~ Xmb ~ 9-218; single Gaussian - 

Xlb = 0-903, Xmb = 4.847; hot-clumpy cuspy - xis = 0-777, Xmb = 4.332; hot-clumpy Gaussian - Xlb = 3.204, Xmb = 3.360. 
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Fig. 3. — Logarithmic i^rms profiles for individual N = 10® GADGET simulations with Qo = l-O- Again, the radial coordinate is scaled by 
the half-mass radius (vertical dashed line). As before, the initial conditions are specified in each panel, and the errorbars indicate the data 
values from the simulations. The thin lines that appear below the errorbars show the behaviors of Urms,r (the lowest line) and the nearly 
identical Urms,0 and Vj-ms,<f)- From these components, one can see the mild tangential velocity anisotropy that exists in these simulated 
equilibria. Two sets of model curves are shown superimposed with the data. As indicated in the legend in the upper-left panel, there are 
isotropic and anisotropic model results. Isotropic model profiles reach nearly constant values near the centers of systems, while anisotropic 
model profiles provide better representations of the data for smaller r. In general, anisotropic LB models produce the smallest velocity 
x? values: single cuspy - Xlb.Iso = 1-108- XLB.aniao = 0-036, Xmb.Iso = 1-072. XMB,aniso = 0-050; single Gaussian - XLB,iso = 1-087, 
XLB.aniso = 0-078, x|,B,isc = 1-778, x|lB,anise = 0-832; hot-clumpy cuspy - X?,B,iso = 0-466, xL,anise = 0-016, xilB.ise = 0-500, 
Xmb, anise = 0-087; hot-clumpy Gaussian - XuB.ise = 0-447, Xlb, anise = 0-038, XMB.ise = 0-325, Xmb, anise = 0-033. Goupled with the 
generally better density fits provided by the LB models, these velocity fits suggest that the LB models are the superior description of the 
phase-space distributions of these simulations. 
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4.— Logarithmic Vrms profiles for individual N = 10^ GADGET simulations with Qq = 0.7. The panels here are analogous to 
Figure [3] The ineffectiveness of isotrc^ic models remains evident, and the superiority of anisotropic LB models is even clearer 
the Qo = 1.0 cases: single cuspy - Xlb.Iso = 2-231, XnB.aniso = 0-214, Xmb.Iso = 1-024, XMB.aniso = 0-644; single Gaussian 
:o = 0-953, xL, anise = 0-072, X^lB.iso = 1-393, XMB.aniso = 1-^24; hot-clumpy cuspy - X?,B,iso = 0-696, X?,B,aniso = 0-043, 
= 1.866, xiiB.aniso = 3-254; hot-clumpy Gaussian - x^B.iso = 0-324, XMB.aniso = 0-167, x|iB,iso = 0-546, xliB.aniso = 1-130- 
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Fig. 5.— Logarithmic density profiles for averaged N = 10^ GADGET simulations with Qq = 1.0. The panels are analogous to those 
in Figure^ The errorbars marking the data points are smaller and the data points show less point-to-point variation than in Figure fTI 
As before, the LB models tend to provide better descriptions of the data: single cuspy — Xlb ~ 4.569, Xmb “ 4.438; single Gaussian - 
Xlb “ 0.635, Xmb “ 3.182; hot-clumpy cuspy - Xlb “ 0.769, Xmb “ 0.652; hot-clumpy Gaussian - Xlb “ 1-389, Xmb “ 0.754. 
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Fig. 6.— Logarithmic density profiles for averaged N = 10^ GADGET simulations with Qq = 0.7. The panels are analogous to those 
in Figure [2] and the data profiles present smoother versions of the outer-profile concavity features seen there. Again, the LB models 
tend to describe the data behavior more completely than the MB models: single cuspy - ~ 7.658, Xmb ~ 44.058; single Gaussian - 

Xlb “ 3-346, Xmb “ 28.463; hot-clumpy cuspy - Xlb “ 0.788, Xmb “ 7.670; hot-clumpy Gaussian - Xlb “ 2.215, Xmb “ 9-064. 
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Fig. 7.— Logarithmic Vrms profiles for averaged N = 10^ GADGET simulations with Qq = 1.0. The panels are analogous to those in 
Figure [ 3 ] These smoother versions again indicate that models incorporating the velocity anisotropy present in the simulations are better 
suited to describing the data. Of the models considered here, the anisotropic LB models provide the best representation of the data: single 
cuspy - XLB.iso = 5-866, XLB.aniso = 0-3™, XMB.iso = 5-629, XMB.aniso = 0-424; single Gaussian - xL.iso = 6-326, Xm, anise = 0-121, 
XMB.iso = 9-294, XMB.aniso = 3-013; hot-clumpy cuspy - XuB.iso = 2-521, XlB, anise = 0-080, XMB.ise = 3-319, Xmb, anise = 0-567; 
hot-clumpy Gaussian - xL.iso = 3-137, Xlb, anise = 0-087, Xmb, iso = 3-531, XMB.aniso = 0-357. 
















Entropy Production 


15 



Fig. 8 .— Logarithmic tirms profiles for averaged N = 10® GADGET simulations with Qq = 0.7. These panels are analogous to 
those In Figure The smaller errorbars in these simulations result In some poorer fits between the anisotropic LB model and the 
data {e.g.,m the clumpy Gaussian simulation, lower-right panel): single cuspy - Xlb iso “ 8.069, Xlb aniso “ 0.730, Xmb iso “ 3.699, 


^MB.a 

X^B.iso = 2.779 


XlvIB.i! 


= 1.029; single Gaussian - X^B.iso = '^•240, Xlb, a 
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Fig. 9.— Logarithmic density profiles for individual N = 10® GADGET simulations with Qq = 1.0. These panels are analogous to those 
in Figure^ The point-to-point variations seen in the single simulations can be relatively large, leading us to question the appropriateness 
of our estimated data point uncertainties. However, the clumpy simulation results are quite smooth, and at least on a relative basis, the 
LB models describe the data better than the MB models: single cuspy - Xlb ~ 15.891, Xmb ~ 16.596; single Gaussian - Xlb “ 6.956, 
^MB “ 7.218; hot-clumpy cuspy - Xlb “ 0.844, Xmb “ 0.859; hot-clumpy Gaussian - Xlb “ 0.561, Xmb “ 1-364. 
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Fig. 10.— Lc^arithmic density profiles for individual N = 10® GADGET simulations with Qq = 0.7. These panels are analogous to 
those in Figure [2] Again, the data points for clumpy simulations show a smoothness not present in the single simulations. Except for the 
single cuspy simulation, the LB models clearly provide a superior description of the data: single cuspy - x^b “ 81.184, Xmb ~ 105.158; 
single Gaussian - Xlb “ 3.572, Xmb “ 32.259; hot-clumpy cuspy - Xlb ~ 1-999, Xmb ~ 28.178; hot-clumpy Gaussian — Xlb ~ 2.730, 
xliB = 22.357. 
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Fig. 11. — Logarithmic Urms profiles for individual N = 10® GADGET simulations with Qo = 1-0. These panels are analogous to those 
in Figure [ 3 ] The models with velocity isotropy continue to be poor descriptors of the data, while the anisotropic LB models provide the 
best fits to the simulation results: single cuspy - Xlb.Iso = 16.044, XLB.aniso = 2-962, Xmb.Iso = 15-723, XMB.aniso = 3-188; single Gaussian 
- xL.isc = 12-313, XLB.aniso = 3-004, xiiB.iso = 12-135, xliB,anise = 3-733; hot-clumpy cuspy - = ^-312, X?,B,aniso = 0-205, 

XMB.iso = 6-728, xiiB.aniso = 0-^27; hot-clumpy Gaussian - x^e.iso = 3-574, XLB.aniso = 0-153, X^ib.Iso = ^-383, xliB.aniso = 1-226- 
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Fig. 12.— Logarithmic i^rms profiles for individual N = 10® GADGET simulations with Qo = 0-7. These panels are analogous to those 
in Figure [4| Again, the anisotropic LB models provide the best descriptions of the data (even though it appears to under-predict the 
single cuspy simulation results for larger r): single cuspy - XLB.iso = 12.735, XLB.aniso = 1-824, XMB.iso = 9-868, XMB.aniso = 2-985; 
single Gaussian - Xlb.Iso = 17-030, XLB,aniso = O-l’^O, Xmb.Iso = 7-721, XMB.aniso = 2.673; hot-clumpy cuspy - Xlb.Iso = 13-464, 
XLB.aniso = 0-806, XMB.iso = 5-429, XMB.aniso = 2-756; hot-clumpy Gaussian - XLB.iso = 15-484, XLB.aniso = 0-756, XMB.iso = 7.722, 
XMB.aniso = 2.308. 
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Fig. 13.— Logarithmic density profiles for individual N = 2^^ GPU-NBODY-6 simulations with Qo = 1.0. These panels are analogous to 
those in Figure^ Evolving systems with this non-softened code produces equilibria with density profiles that are nearly indistinguishable 
from those discussed earlier. As with the GADGET results, LB models provide better representations of the data: single cuspy - Xlb “ 
1.952, Xmb “ 2.064; single Gaussian - Xlb “ 0.390, Xmb “ 1-537; hot-clumpy cuspy - Xlb “ 1-096, Xmb “ 0.964; hot-clumpy Gaussian 
- xIb = 0-595. xIiB = 0.584. 














Entropy Production 


21 



Fig. 14.— Logarithmic density profiles for individual N = 2^^ GPU-NBODY-6 simulations with Qq = 0.7. These panels are analogous 
to those in Figure [2] The fact that the LB models can reproduce the slight concave features in the outer regions of the profiles again lead 
to them being preferred to the MB models: single cuspy - x^b ~ 3.653, Xmb “ 17.118; single Gaussian - Xlb ~ 1-420, Xmb ~ 13.025; 
hot-clumpy cuspy - Xlb “ 3.191, Xmb “ 13.159; hot-clumpy Gaussian - Xlb “ 0-735, Xmb “ 8-175. 
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Fig. 15.— Logarithmic t;rms profiles for individual N = 2^^ GPU-NBODY-6 simulations with Qq = 1.0. These panels are analogous 
to those in Figure [3] The trend for anisotropic LB models to provide the best description of the Urms data continues. Note that the 
clumpy cuspy simulation results in a nearly isotrcy^ic velocity distribution, and that the inner part of the profile is decently well-described 
by the isotropic LB prediction: single cuspy - xtB.iso = 1-741, XLB.aniso = 0-048, XMB.iso = 1-639, XMB.aniso = 0-058; single Gaussian 
- XLB.iso = 0.904, XLB.aniso = 0-032, Xmb.Iso = 1-^34, XMB.aniso = 0-953; hot-clumpy cuspy - XuB.iso = 0.508, XuB.aniso = 0.194, 
XMB.iso = 0-711, XMB.aniso = 0-310; hot-clumpy Gaussian - XuB.iso = 1-278, Xlb, anise = 0-192, Xmb.Iso = 1-787, Xmb, anise = 0-560. 
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Fig. 16.— Logarithmic t;rms profiles for individual N = 2^^ GPU-NBODY-6 simulations with Qq = 0.7. These panels are analogous 
to those in Figure [4] The anisotropic LB models continue to represent the data behavior better than the other models considered: single 
cuspy - XLB.iso = 1-899. XLB.aniso = 0-128, XMB.iso = 0-863, XMB.aniso = 0-725; single Gaussian - Xlb.Iso = 1-956, Xlb, anise = 0-114, 
XMB.iso = 1-365, xllB,anise = l-^^S; hot-clumpy cuspy - X?,B,ise = 0-516, X?,b, anise = 0-147, xllB.ise = 0-665, xllB,anise = 1-177; 
hot-clumpy Gaussian - x^B.ise = 1-560, X?,b, anise = 0-384, x^B.ise = 0-877, X^ib, anise = 1-306- 



is calculated from the phase-space macro-cell occupation values (Equation [Jll, with a value of i^lb = 10^. Two different macro-cell volumes 
(indicated in the legend) have been used to create comparison curves. Adopting a smaller macro-cell volume increases the zero point of 
the curve and reduces the small-scale variations in Slb but leaves the overall behavior unchanged, (b) The virial ratio Q as a function of 
time for the same simulation. The initial increase in S'lb appears to coincide with the initial growth in Q. 
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Fig. 18.— The macroscopic entropy production rate versus time in the individual single, cuspy N = 10® GADGET simulation with 
Qq = 0.7. The production rate is calculated by determining gradients in the kinetic temperature and mean velocity field (Equation [8jl. 
The rather coarse grid used to determine the gradients may play an important role in the essentially null result shown. If the macroscopic 
entropy production could be properly calculated, one would expect a rather large positive spike in the interval 0 < t/T < 5. 



















